#!/usr/bin/python
"""This module calculates correlation coefficient.

This module takes output file, containing predicted sites and calculates correlation
coefficient with real frequency matrix.
"""

# by Vineet Gupta <vineet@ural.wustl.edu>
# cal_cor.py 2-2-2007 15:05

# TO DO:
# 

import sys
from alignment import Alignment

def pearson(X, Y):
    '''Calculates Pearson correlation between elements in list X & Y.'''

    if len(X) != len(Y):
        print "Attemp to calculate correlation between dissimilar number of counts."
        sys.exit(1)

    N = len(X) # N is total number of elements in list X & Y
    sigma_x = 0.0
    sigma_y = 0.0
    sigma_x2 = 0.0
    sigma_y2 = 0.0
    sigma_xy = 0.0

    for i in range(0, N):
        sigma_x += X[i]
        sigma_y += Y[i]
        sigma_x2 += X[i]*X[i]
        sigma_y2 += Y[i]*Y[i]
        sigma_xy += X[i] * Y[i]

    num = (sigma_xy) - ((float)(sigma_x) * sigma_y / N)
    deno = pow((sigma_x2 - pow(sigma_x, 2) / (float)(N)) * (sigma_y2 - pow(sigma_y, 2) / (float)(N)), 0.5)

    return (num / deno)

def main():
    a = Alignment(sys.argv[1])
    # True matrix for half site binding for Mnt protein
    hs =  [0.14, 0.06, 0.15, 0.23, 0.40, 0.10, 0.21, \
    0.07, 0.10, 0.01, 0.08, 0.30, 0.79, 0.43, \
    0.75, 0.14, 0.76, 0.63, 0.10, 0.03, 0.15, \
    0.04, 0.70, 0.08, 0.06, 0.20, 0.08, 0.21]
    
    print pearson(a.cor_matrix(), hs)
    
if __name__ == '__main__':
    main()
    